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This  paper  deals  with  the  development  of  a  mathematical  foundation  and  practical  algorithms  to 
compute  the  safety  fans  to  he  used  in  conjunction  with  aerodynamic  deceleration  payload  delivery 
systems  tests.  The  paper  provides  formulas  to  compute  a  descent  /  gliding  trajectory  for  a 
multiple-stage  system  and  then  considers  several  cases  of  major  failures.  It  also  introduces  the  concept 
of  hallistic  winds,  simplifying  calculation  of  the  safety  fans  so  that  the  computer  is  no  longer  needed, 
and  shows  how  to  carefully  calculate  these  hallistic  winds  based  on  the  best  known  winds  aloft.  Finally, 
the  paper  introduces  the  graphical  user  interface  developed  and  tested  by  the  authors  that  is  supposed 
to  provide  test  planning  officers  with  a  unique  variety  of  trade-off  options  as  well  as  to  effectively 
represent  all  safety-related  information  both  to  the  ground  personal  and  to  the  crew  performing  a 
payload  delivery  from  a  transport  aircraft  or  helicopter.  The  paper  ends  with  conclusions. 
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Notations 

=  altitude,  speed  and  direetion  of  the  aireraft  at  the  release  point 
=  altitude  AGL,  at  whieh  the  2'“*  and  3'^'*  stages  are  deployed  (if  present) 

=  altitude  loss  before  the  1®‘  stage  canopy  deployment 
=  system’s  Northing,  Easting  and  altitude 

=  Northing  and  Easting  at  the  altitude  of  the  /-th  stage  deployment  (1=1,2, 3) 
=  speed  and  direction  of  the  wind  at  altitude  h 
=  descent  rate  and  glide  ratio  during  the  1-th  stage  (1=1, 2, 3) 


1.  Introduction 

AS  with  any  other  technical  system,  aerodynamic  delivery  systems  (ADS)  are  subject  to  failures.  That  is  why 
numerous  tests  are  required  to  tune  the  system,  make  it  more  reliable  and  therefore  minimize  the  risk  of  failure. 
A  variety  of  ADSs  rely  on  a  parachute,  a  cluster  of  parachutes,  or  most  recently  autonomous  parafoils  to  deliver 
payloads  of  different  weight  m  (up  to  several  tons)  from  different  altitudes  (up  to  several  kilometers).  The 

canopy  opening  failures  result  in  a  ballistic  trajectory  that  transforms  the  initial  kinetic  (O.SmVf^)  and  potential 
( ingHg )  energy  into  kinetic  energy  at  impact. 

A  failed,  ballistically  falling  system  can  cause  huge  damage  at  impact.  That  is  why  it  is  important  to  foresee  such 
an  unfortunate  event  and  take  certain  precautions  organizing  these  drops  in  remote  and  uninhabited  areas.  Basic 
physics  provides  a  very  simple  formula  to  calculate  the  ground  distance  the  failed  ADS  can  cover.  If  we  neglect 
aerodynamic  drag  (yielding  the  largest  ground  distance  estimate  possible)  this  formula  is  as  simple  as  V^T ,  where 

the  time  of  fall  T  can  be  estimated  as  T  =  ^j2Hgg~'  .  Of  course  if  we  have  strong  horizontal  winds,  they  will 

definitely  affect  this  range  (horizontal  projection  of  the  trajectory  will  in  this  case  not  be  a  straight  line  any  longer), 
though  the  effect  will  be  small  as  the  ADS  would  fall  through  these  winds  fairly  quickly. 

In  case  of  partial  canopy  opening  failures,  however,  the  effect  of  winds  will  be  stronger.  This  is  not  only  because 
the  system  will  descend  slower  and  take  longer  to  reach  the  ground,  but  also  because  we’ll  have  to  account  for 
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aerodynamic  drag.  In  the  case  where  we  are  dealing  with  a  failed  parafoil  the  situation  can  be  even  worse  because  of 
the  high  glide  ratio  of  such  systems.  With  a  fully  deployed  canopy  but  with  the  controls  stuck  in  some 
predetermined  position  or  in  the  case  of  a  controls  failure,  the  parafoil-based  ADS  can  glide  far  away  from  the 
intended  point  of  impact  (IPI). 

This  paper  develops  a  mathematical  foundation  and  corresponding  software  for  reliable  computation  of  safety 
fans  for  multistage  systems  (when  the  deployment  of  multiple  canopies  is  staged)  and  is  organized  as  follows. 
Section  II  provides  the  mathematical  foundation  for  computing  safety  fans  for  multistage  systems  in  the  different 
types  of  failure  modes  described  above  (canopy  opening  failure  and  controls  failure  for  controllable  ADSs).  Next, 
Section  III  presents  a  detailed  example  of  computing  these  safety  fans  for  a  two-stage  system.  Section  IV  introduces 
the  concept  of  so-called  ballistic  winds,  the  instrument  allowing  major  simplification  of  safety  fan  computation.  This 
simplification  provides  a  reliable  method  for  a  range  safety  officer  to  compute  these  fans  even  without  a  computer 
by  using  a  simple  look-up  table  prepared  upfront  using  the  best  known  winds  at  the  drop  zone  (DZ).  Section  V 
presents  a  graphical  user  interface  (GUI)  developed  in  the  MATLAB  development  environment  that  provides  much 
more  flexibility  in  computing  safety  fans,  suggesting  release  points  and  directions,  and  therefore  allows  for  assessing 
and  mitigating  risks.  The  paper  ends  with  conclusions. 

II.  Computation  of  Safety  Fans 

This  section  considers  two  major  failures.  The  first  case  is  when  the  very  first  stage  canopy  (or  a  single  canopy 
in  a  one-stage  system)  doesn’t  open  at  all  and  the  ADS  falls  ballistically.  This  will  be  the  canopy  opening  failure 
case.  The  second  case  is  when  the  first  stage  is  deployed  successfully  but  there  are  certain  problems  with  steering  the 
system  in  the  desired  direction.  This  case  will  be  referred  to  as  the  controls  failure  case.  The  effect  of  winds  in  both 
cases  will  be  accounted  for  as  described  above. 


A.  Canopy  Opening  Failnre 

Assuming  that  aerodynamic  drag  can  be  neglected,  the  equations  driving  the  ADS  motion  are: 

dx 

—  =  w{h)cosy/^{h)  +  Vf^cos'¥ 
dt 


dy 

dt 


=  w{h) sin y/„{h)  +  sin 'Pj 


(1) 


dh 

-r  =  8t 
dt 


Here  the  vertical  motion  is  driven  by  the  gravity  force,  while  the  horizontal  motion  is  caused  by  the  initial  velocity 
and  variable  winds.  Equations  (1)  can  be  further  reduced  to: 


dx  dh 

- =  w{n) cosiy ^{n)  +  cos'V 

dh  dt 

dy  dh 

dh  dt 

dh 

dt 


(2) 


From  the  last  equation  it  follows  that 


(3) 


With  this  relationship  between  elapsed  time  and  altitude  the  first  two  equations  in  (2)  can  be  further  rewritten  as 

dx  _  w(/?)cost//„(/?)-i-FoCos'Fo 
dh~  pgiH,-h) 

dy  w(/?)sin(//^^(/!)-i-F(,sin'F„ 
dh~  ^2g{H,-h) 

Now  the  X  andy  coordinates  at  any  arbitrary  altitude  h  can  be  obtained  via  integrating  Eqs.(4): 


(4) 
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(5) 


cos (//„(/?) +  F(, cos 'Po 

x{n)  =  Xg  +  - 1  = - an 

i  pgiH,-h) 

. . ,  ,  ?  sin v^„(h)  +  Vo  sin '¥ „ 

yW  =  Po  +  - 1  ,  - dh 

i  pgiHo-h) 

The  integrals  in  these  equations  can  be  further  split  into  two  separate  components,  so  that 


{  j2g(N,-/,)  i 

,  ?  w(/?)  sin  (//„(/;)  [ 

yW  =  To  +  I  I  ,  ydh  +  . 

i  pgiHo-h)  V 


H  -h 

2—2 — -F„cos'P„ 


H  -h 

2^— FgSinT'g 

g 


(6) 


If  the  parameters  w{h)  and  are  given  as  a  look-up  table  with  the  elements:  and  k  =  \,...,N , 

(K  >  )’  then  the  effect  of  the  winds  can  be  reduced  to  the  following  finite  sums: 


1 


•  w(h)  cost// ^(h) 


tn  pg{Ho-h,) 


1 


•  w{h)siny/^{h) 


dh^Y 


(7) 


pg{Ho-h)  pg(Ho-h,) 


(where  //„  corresponds  to  and  h  to  h^  ). 

Next,  let  us  consider  the  case  when  the  canopy  opens  successfully  but  the  controls  are  stuck  in  some  fixed 
position,  causing  the  ADS  to  glide  with  a  constant  heading  (ADS  dynamics  can  be  neglected). 

B.  Controls  Failure 

A  multiple-stage  ADS,  say  with  three  stages,  assumes  the  following  consecutive  descent  phases: 

•  Ballistic  trajectory  before  the  first  stage  is  fully  deployed; 

•  Deployment  of  Stage  1 ; 

•  D  ep  loyment  o  f  Stage  2 ; 

•  Deployment  of  Stage  3 . 

These  stages  are  described  below. 

Ballistic  trajectory  before  the  first  stage  is  fully  deployed.  We  may  apply  Eqs.(6),  which  for  the  altitude  range 
of  hG[Ho-AH;Ho]  yield 


x{h)  =  Vg  -I-  I 


•  w(h)  cos  t//^(h) 


pgiHo-h) 


dh  +  2 


H„-h. 


g 


-FgCOST'g 


ng  ,  r  w{h)smii/„{h) 

y{h)  =  yo+  I  ,  , 

H„-Mf  ^2g[Ho  -h) 


dh  +,  2 


H.-h 


(8) 


g 


VoSm'¥o 


Specifically,  by  the  time  (altitude)  the  first  stage  deploys  we  have 


1 


•  w(h)  cos  t//^(h) 


pgiH.-h) 


dh+\2 


AH. 


g 


-F„  COSTCO 


Ti  =  To  +  j 


^  w{h)smy/^{h) 


pgi^o-h) 


dh  +,  2 


AH 


(9) 


g 


FgSinT'g 


From  this  point  on  we  assume  that  the  initial  horizontal  speed  caused  by  deploying  from  an  aircraft  flying  with 
speed  Vq  is  phased  out  completely.  For  simplicity  we  will  not  consider  the  gradual  transition  from  a  ballistic  fall  to 
the  normal  gliding  descent.  Indeed,  accounting  for  this  transition  may  cause  some  minor  changes  in  computation  of 
the  safety  fan  but  being  on  the  conservative  side,  we  will  not  care  that  the  size  of  the  fan  might  be  smaller  without 
the  transition.  In  any  case,  it  is  likely  that  wind  changes  may  introduce  lager  error  than  neglect  of  the  transition. 
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Descent  of  a  fully  deployed  system.  After  canopy  deployment  the  trajectory  obeys  the  following  set  of  ordinary 
differential  equations: 

dx 

—  =  w(h)co?,y/^(h)  +  v^fiR^  cost// 


dt 

—  =  w{h)  sin  (//„(/?)  +  v^,G7?j  sin  y/ 
dt 

dh 


(10) 


Here  we  assume  that  the  descent  rate  does  not  depend  on  the  altitude  (air  density),  and  is  constant  during  the  entire 
decent.  Using  the  third  equation  as  demonstrated  previously,  the  first  two  equations  can  be  reduced  to 

dx 

—  -  cos  ¥w(h)  +  ^^1  cos  y/ 


dh 

—  =  v~lw{h) sin (//„(/?)  +  GR^  sin y/ 
dh 

Now,  we  can  integrate  these  equations  with  respect  to  an  altitude,  so  that 

x{h)  =  x^+  j  {v~lw{h)cosy/^{h)  +  GR^cosy/^dh 

h 

H„-AH 

y{h)-y^+  j  {v~lw{h)^,my/^{h)  +  GR^^,my/^dh 

h 

Again,  the  integrals  can  be  split  into  two  parts 

H„-AH 

x{h)  =  Xi+  j  v^' w(h)  cos  y/„(h)dh  +  (Hg- AH -h)GR^  cosy/ 


(11) 


(12) 


h 


(13) 


y{h)  =  y^+  j  v2'w{h)smy/^{h)dh  +  {H^- AH -h)GR^smy/ 


where  the  third  terms  represent  the  contribution  from  the  stuck  controls  and  the  second  terms  provide  the 
contribution  of  the  variable  winds.  By  the  altitude  of  deployment  of  the  second  stage  we  will  have 


Vj  =  Vj  +  j  v~^w(h)  cos  y/^  (h)dh  +  -  AH  -  )GR^  cos  y/ 


Hi 
H„-AH 


(14) 


j  v^l'w{h)smy/^{h)dh  +  {Hf^- Mi  -  H^)GR^smy/ 


Deployment  of  Stage  2.  If  Stage  2  is  present  then  the  equations  of  motion  will  be  similar  to  those  of 
Eqs.(12)-(14).  Specifically,  we  can  write 

Hi 

x(h)  =  X2+  [  v'2v/{h) cosy/ ^{h)dh  +  {H^-h)GR^  cosy/ 

I  (15) 

>’(^)  =  >’2  +  j  v~lw{h)smy/^{h)dh  +  {H^-h)GRjSmy/ 


and 


“2 

Xj=X2  +  j  v^^w(h) cos y/^(h)dh  +  {H^  - H, )G7?2  cos y/ 

Hi 

Hi 

y3^y2+  j  Ui‘ w(h)smy/^(h)dh  +  (H^  - H^ ) GR^  sin y/ 


(16) 


respectively. 

Deployment  of  Stage  3.  Finally,  for  Stage  3  (if  present)  we  will  have 
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(17) 


x(h)  =  x^  +  j  v~lw{h)co^y/^{h)dh  +  {H^  -h)GR^  cosy/ 

h 

y(h)  =  +  j  v;lw{h) smy/^{h)dh  +  {H^  - h)GR^  sin y/ 

h 

and 

x(0)  =  X3  +  j  v~2w{h)cosy/^,Qi)dh  + HfiR^  cosy/ 

0 

j(0)  =  J3  +  j  v~^w(h)smy/^{h)dh  + H^GR^  sin(// 

0 

Now  that  we  have  developed  the  equations  of  motion  for  a  multi-stage  system,  let  us  eonsider  a  representative 
example. 

Ill,  Example  of  Safety  Fans  for  a  Two-Stage  System 

Consider  a  system  that  has  two  stages,  whieh  are  deseribed  as  follows: 

Vi,  mis  GRi  /Url  {H^),m 

Stage  1  18  2.5  666 

Stage  2  6  0  333 

Assume  we  release  it  from  =  1, 500m  at  -  60m  /  s  due  North-West.  The  winds  are  given  in  a  look-up 
table  and  stored  in  the  Winds  ( : ,  1 :  3 )  matrix  (where  the  first  eolumn  eontains  altitude,  seeond  -  wind  direetion, 
and  third  -  wind  magnitude).  (Hereafter  we  will  use  fragments  of  self-explanatory  aetual  eode  written  in 
MATLAB.) 

Depending  on  the  data  eolleetion  method  the  wind  measurements  may  be  not  distributed  evenly  in  altitude.  So, 
before  eomputing  definite  integrals  in  Eqs.(6),(8),(12),(15),  or  (17)  numerieally  (using  trapezoidal  numerieal 
integration)  the  original  look-up  table  /i^ ,  ,  y/^^,/.  should  be  approximated  with,  for  example,  pieeewise  eubie 

Hemiite  interpolating  polynomials  w{h)  and  y/„{h)  : 

%%  Setting  initial  conditions 

V0=60;  H0=1500;  Psi0=45*pi/180;  g=9.81;  DR=[18;  6];  GR=[2.5  0]; 

DH=666;  H2=333; 
h=linspace (0, HO, 98) ; 

h (end+1 ) =H0-DH;  h (end+1 ) =H2 ;  h=sort  (h) ; 

%%  Producing  spline  interpolations  for  the  winds 
psih  =  pchip (Winds (:, 1 ), Winds (:, 2 ), h) *pi/180 ; 
wh  =  pchip (Winds ( ; , 1 ) , Winds ( ; , 3 ) , h) ; 

Next,  the  ballistie  trajeetory  due  to  initial  veloeity  (the  last  terms  in  Eqs.(6))  ean  be  eomputed  as  follows: 

%%  Computing  a  ballistic  fall  if  no  winds  were  present 
xtl=sqrt (2* (HO-h) /g) *V0*cos (PsiO) ; 
ytl=sqrt (2* (HO-h) /g) *V0*sin (PsiO) ; 

If  aeeounting  for  the  winds  (the  seeond  terms  in  Eqs.(6))  the  trajeetory  beeomes 

%%  Computing  contribution  of  winds 
for  i=l : 9  9 

xw(i)  =  trapz(h(i:100),wh(i:100).*cos(psih(i:100)). /sqrt (2*g*H0-h(i:100) )  )  ; 
yw(i)  =  trapz(h(i:100) ,wh(i:100) .*sin(psih(i:100) ) . /sqrt (2*g*H0-h(i:100) ) ) ; 
end 

xt2=xtl+[xw  0]; 
yt2=ytl+ [yw  0 ] ; 

Both  trajectories  are  shown  in  Fig.l  (the  star  represents  a  release  point).  This  latter  trajeetory  represents  the  ballistie 
fall  if  the  first-  (and  second-)  stage  canopy  does  not  open.  The  winds  are  pushing  the  trajectory  to  the  South-East 
but  because  the  system  falls  ballistically  (very  quickly)  the  winds  introduce  no  major  change. 
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Fig.l  Ballistic  trajectories  with  no  wind  (a)  and  winds  (b). 

Now,  deploying  the  first  stage  canopy  means  that  we  are  slowing  down  the  descent  and  the  winds  will  produce  a 
much  greater  effect.  Fig. 2a  shows  the  addition  of  the  second  terms  in  Eq.(14)  computed  as  follows: 

%%  Deploying  the  1st  stage 
indexl=f ind (h>H0-DH) ;  Indl=indexl ( 1 ) -1 ; 
for  1=1 : Indl-1 

xSl(i)  =  trapz(h(i:Indl), wh (i;Indl)  .*cos  (psih(i:Indl) )  ./DR(1) ) ; 
ySl (i )  =  trapz  (h  (i : Indl ) , wh (i;Indl)  .*sin(psih(i:Indl) )  ./DR(1) ) ; 

end 

xSl  =  xSl+xt2 ( Indl ) ;  ySl  =  ySl+yt2 ( Indl ) ; 

The  only  difference  from  Eq.(14)  is  that  the  trajectory  is  propagated  all  way  down  to  the  surface  rather  than  to  H2. 
When  the  second  stage  is  deployed  (which  leads  to  even  slower  rate  of  descent)  the  winds  play  an  even  greater  role 
(Fig.2b): 

%%  Deploying  the  2nd  stage 
index2=f ind (h>H2 ) ;  Ind2=index2 ( 1 ) -1 ; 
for  1=1 : Ind2-1 

xS2(i)  =  trapz(h(i: Ind2 ) , wh (i;Ind2) .*cos(psih(i:Ind2)) ./DR (2)) ; 
yS2(i)  =  trapz(h(i: Ind2 ) , wh (i:Ind2)  .*sin(psih(i:Ind2))  ./DR  (2)) ; 
end 

xS2  =  xS2+xSl  (Ind2)  ;  yS2  =  yS2  +  ySl ( Ind2 ) ; 


Fig.2  The  first-stage  opening  (no  glide)  (a)  and  the  seeond-stage  opening  (no  glide)  (b). 

Finally,  the  effect  of  the  controls  failure  mode  (controls  stuck  in  one  position  leading  to  a  glide  with  constant 
ground  track  angle)  is  reduced  to  drawing  a  circle  of  appropriate  radius  (the  last  term  in  Eq.(14))  around  the  last 
point  of  the  current  stage. 
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%%  Computing  a  safety  fan  for  the  first  stage 
[X,Y,Z]  =  cylinder ( (H0-DH-H2) *GR(1) ) ; 

%%  Computing  a  safety  fan  for  the  second  stage 
[X,Y,Z]  =  cylinder ( (H0-DH-H2) *GR(1) +H2*GR (2) ) ; 

The  results  are  presented  in  Fig. 3.  Sinee  the  speeifie  system  we  are  eonsidering  uses  round  paraehutes  with  no  glide 
as  a  seeond  stage,  the  final  safety  fan  is  simply  first-stage  safety  fan  shifted  South-East. 


I _ I _ I _ I _ I,..  „  I _ I _ I _ 

-1000  -500  0  500  1000  1500  2000  2500  3000 


Easting,  m 

Fig.3  Safety  fans  in  case  the  controls  are  stock  in  a  fixed  position. 

Computation  of  the  safety  fans  for  different  multistage  systems  is  earried  out  similarly. 

As  seen,  the  eomputations  above  are  not  very  eomplieated  and  yet,  they  are  unique  for  eaeh  system  and  earried 
out  numerieally  based  on  the  unevenly-spaeed  measurements  of  the  winds  aloft.  If  two  different  systems  are  to  be 
dropped,  these  ealeulations  should  be  performed  for  eaeh  system  individually.  The  goal,  however,  is  to  have  another 
teehnique  allowing  the  use  of  some  sort  of  normalized  winds  that  ean  be  eomputed  upfront  onee  and  then  be  applied 
to  any  system.  This  brings  us  to  the  next  seetion,  whieh  diseusses  the  so-ealled  ballistie  winds. 


IV.  Computation  of  Ballistic  Winds 


For  operational  use  the  vertieal  winds  profile  ean  be  redueed  to  the  so-ealled  ballistie  winds.  In  other  words,  if  at 
some  altitude  H  we  have  a  ballistie  wind  of  magnitude  W  and  direetion  'i'w,  then  the  effeet  of  variable  winds  for 
some  system  with  deseent  rate  Vy  on  its  way  from  altitude  //down  to  the  surfaee  is  redueed  to  simple  formulas: 


x{h)  -  — IF  cos 
y(h)  =  yWsm^‘^ 


(19) 


A.  Deriving  Ballistic  Winds 

Comparing  Eqs.(19)  with  Eqs.(18)  it  becomes  clear  that 

?  1  H 

J  — w{h)c,osy/^{h)dh  =  — W  cos^F^ 


.Vy 


^  1  H 

I  —  w{h)smy/^{h)dh  =  — IF  sin 'F^ 


(20) 


.Vy 


Substituting  integrals  with  the  finite  sum  of  trapezoids  based  on  the  discrete  values  of  and  ,  k  =  I,...,  M  , 

we  get 
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(21) 


^k-\)  ~  « 


k=2 

M 


^A-l  )  ~ 

*=2  ^ 

The  index  starts  from  2  because  by  definition  the  winds  measurements  at  the  lowest  altitude  can  be  considered 
ballistic  winds  at  this  altitude. 

From  Eqs.(21)  it  further  follows  that 


Z(^*  -^A-i)(wt  sin +vi;,_|  smy/„.^) 


tan'F 


W;M  M 


w„  = 


Z(^A  -K-^)[^k  cosy/^..,_,) 

j(  M  7~M  Y 

J  Z(^A“Vi)(ncos(/^„,.,+vi;4_iCOS(/^„,,_,)  +  X(^A“Vi)(nsmt/^„,*+w,_,sin(/^„,.,_,) 

VVA-=2  j  Va=2  J 


A=2 

1 


(22) 


For  the  specific  case  when  hi.  -  hi._i  =  Ah  =  const ,  k  =  2,...,M  ,  Eqs.(22)  can  be  further  reduced  to 


M 

Z(^A  sin +  w,_;  sin 


tan'F^^.„  = 


W-,M  M 


Z(v^A  ^ow„,k  +  Wt-i  cosy/^..,_i ) 


(23) 


Wa  +  w,_i  cos I  +1  X(^A  sin?/^,,*  +  sin ) 


B.  Example  of  Computing  Ballistic  Winds 

Let  us  consider  a  simple  example.  Suppose  that  the  winds  aloft  were  measured  as 


Alt,  ft 

Direction 

0 

93 

1000 

101 

2000 

119 

3000 

129 

4000 

135 

5000 

139 

6000 

145 

7000 

150 

8000 

152 

9000 

151 

10000 

144 

deg  Magnitude,  kts 

10 
36 
25 
21 
19 
17 
16 
16 
16 
16 
16 


and  stored  into  the  Winds  matrix.  Then,  the  following  script 

BWinds ( ; ,  1 ) =Winds ( :  ,  1 ) ;  BWinds (1,2) =Winds (1,2)  ;  BWinds (1,3) =Winds (1,3) ; 
for  n=2 : 1 1 

dh=dif f (Winds(l:n,l)); 

wE=sum (dh . * (Winds (2 : n,  3 )  .  *sin (Winds (2 : n, 2 )  *pi/180)  +  ... 

Winds (1 :n-l, 3) . *sin (Winds (1 :n-l, 2) *pi/180) ))/2; 
wN=sum (dh . * (Winds (2 : n, 3 ) .  *cos (Winds (2 : n, 2 )  *pi/180)+... 

Winds(l:n-1,3) .*cos(Winds(l:n-l,2)*pi/180)) )/2; 

BWinds (n, 2 ) =atan2 (wE, wN) *  18  0/pi ; 

BWinds (n, 3 ) =sqrt (wE^2+wN^2 ) /Winds (n, 1 ) ; 
end 


converts  them  into  the  ballistic  winds: 


o 

o 

93.0 

10.0 

1000.0 

99.3 

23.0 

2000.0 

104.4 

26.5 

3000.0 

110.2 

25.0 
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4000.0 

114.7 

23.5 

5000.0 

118.2 

22.1 

6000.0 

121.2 

21.0 

7000.0 

124.1 

20.1 

8000.0 

126.8 

19.4 

9000.0 

129.1 

18.8 

10000.0 

130.6 

18.5 

Now,  if  for  instance  we  want  to  see  what  would  happen  to  the  second  stage  of  the  system  described  in  Section  III 
if  it  were  deployed  at  3,000ft,  then  instead  of  calculating  integrals  in  Eq.(18)  (which  would  result  in  a  trajectory 
shown  in  Fig. 4)  we  could  simply  apply  Eq.(19): 

X=BWinds (3, 1) /GR(2) *BWinds (3, 3) *cos (BWinds (3, 2) *pi/180) ; 

Y=BWinds (3,1) /GR ( 2 ) * BWinds ( 3, 3 ) *sin (BWinds (3,2)*pi/180); 

This  accurately  computes  the  forecast  point  of  impact  (shown  in  Fig.4  as  a  solid  circle)  in  one  shot. 


Fig.4  Calculation  of  a  predicated  impact  point  based  on  ballistic  winds. 

To  summarize,  the  ballistic  winds  (when  computed  correctly  according  to  Eqs.(22)  and  (23))  provide  an 
excellent  tool  to  rapidly  assess  the  impact  point  in  case  of  canopy  opening  failure  and  descent  of  a  symmetric 
canopy  ADS  (with  zero  glide  ratio).  However,  the  descent  of  a  parafoil  ADS  still  requires  rigorous  computations.  As 
a  result,  the  next  section  presents  a  GUI  developed  to  perform  these  computations  behind  the  scene,  while 
visualizing  the  drop  planning  and  all  safety  related  information  on  one  screen.  This  provides  the  capability  to 
conveniently  change  /  adjust  /  play  with  different  parameters. 

V,  Interactive  GUI 

Figure  5  presents  a  screen  shot  of  the  developed  GUI  which  performs  various  computations  and  manipulations 
(the  factual  data  is  replaced  with  fake  data).  This  GUI  was  developed  in  the  MATLAB  development  environment 
using  the  Open  GUI  Layout  Editor  GUIDE.  The  code  that  runs  this  GUI  contains  over  3000  lines  and  the  GUI  itself 
is  subdivided  into  several  panels,  described  next. 

A.  Winds  Panel 

The  Winds  panel  allows  the  user  to  browse  for  wind  fdes  and  chose  the  most  appropriate  one.  The  wind  source 
fdes  are  sounding,  Windpack,  balloon  and  JAAWIN  winds.  Examples  of  heading  and  the  first  three  lines  of  wind 
data  for  each  of  these  files  are  shown  in  Figs. 6-9.  To  read  each  of  these  files  a  corresponding  function  was  created. 
For  instance,  the  following  set  of  commands 

%%  Reading  heading  information 

[Y,M, D,H,fl,f2,f3,lat,lon]=. . . 

textread (FileName, ' %4c-%2c-%2c_%2c  %s  i,j=(  %f  ,  %f  ),lat,long=(  %f  ,  %f',... 

1 , ' header lines 2 ) ; 

Ldate=[M  '/'  D  '/'  Y] ; ,  Ltime=[H  '  ;  00 '  ]  ; 

file  =  textread (FileName,  '%s',  'delimiter',  '\n'); 

%%  Reading  numerical  data 


9 

American  Institute  of  Aeronautics  and  Astronautics 


N=length ( f ile ) ;  i  =  6;  Winds=[0  0  0 ]  ; 
while  i  <=  N 

[q, w, e, r, t, y, u] =textread (FileName, ' %f  %f  %f  %f  %f  %f  %f 1 headerlines i ) ; 
Winds=vertcat (Winds , [w  y  t] ) ;  %  the  format  is  altitude,  Wdir,  Wspeed 
i=i+l; 
end 

Winds (1 , : ) = [ ] ; 

reads  data  from  the  sounding  fde  (Fig. 6). 


Fig.  5  Safety  fans  GUI. 


37  2 

2008-12-16_00Z, DZIDOO  i,j=(  15.9,  15 . 6 ) , lat , long= ( 54 . 832 , 83 . 104 ) 

SFALT,  SFPRES 
381  968.3 


LINE 

AGL  (m) 

T(C) 

RH(%) 

WSPD (m/s) 

WDD 

P  (mb) 

1 

2.0 

11.76 

50.93 

9.33 

182.93 

968.04 

2 

15.2 

11 . 64 

51.16 

10.97 

182 . 91 

966.53 

3 

56.2 

11.26 

52.30 

13.47 

182.37 

961.80 

Fig.  6  Example  of  a  winds  sonnding  file. 


Time (UTC) 

Lat (deg) 

Long (deg) 

Altitude (M-HAE) 

1 

Q  StdDev 

VE (M/sec) 

VN (M/Sec) 

VUp (M/Sec 

53728.000 

54 . 832114 

83.104813 

5613.796 

4 

0.697 

-0.970 

91.957 

-0.030 

53728.100 

54 . 832114 

83.104813 

5613.817 

4 

0.697 

-0.967 

91 .811 

0.142 

53728.200 

54 . 832114 

83.104813 

5613.828 

4 

0.697 

-1.097 

91.665 

0.102 

Fig.  7  Example  of  a  Windpack  winds  file. 
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0700L  FLIGHT 


Station  Name 

:  DZOO 

Radiosonde  Type: 

RS92-SGP 

Latitude : 

54.83 

Radiosonde  Number 

:  DXXXXXX 

Longitude : 

83.10 

Altitude 

:  1322  feet  (  403  Meters) 

Ground  Check 

Corrections 

Pressure:  0. 

44  mb 

RHl;  1.07  % 

Temperature : 

-0.84  °C 

RH2:  -0. 

41  % 

Launch  Time: 

14:14  UTC 

Sounding 

Started  ^ 

on;  12/15/2008 

Dew  Pt 

Air 

Wind 

Wind 

Altitude 

Press 

Temp  RH 

Temp 

Density 

Direct 

Speed 

(FT  above  msl)  (MB) 

(  °F)  (%) 

(°F) 

(g/m^3) 

(degs) 

(Kts) 

1322 

968.0 

45.6  44 

24 . 9 

1197 . 94 

124.0 

9.0 

1400 

965.2 

45.4  41 

23.0 

1195.07 

127.0 

8.9 

1500 

961.6 

45.5  41 

23.3 

1190.34 

129.5 

10.7 

Fig.  8  Example  of  a  balloon  winds  file. 

A.  Forecast 

Time,  Drop 

Zone  Location, 

Elevation 

Valid  Time:  1400Z  27  Oct  2008  Data  Source:  MM5 

Forecast  Model  Grid  Pt:  Latitude->  54.83N  Longitude->  83.104E 

Forecast  Model  Grid  Pt  Elevation:  1719  feet  above  MSL 


Point  of  Impact:  Latitude->  54.83N  Longitude->  83.104W 

Forecast  Model  Grid  Point  Bearing/Range  from  Point  of  Impact:  90DEG/0.00KM 
(*Note!  Model  Grid  Pt  is  nearest  to  entered  Lat/Lon  and  is  in  decimal  degrees  format.) 
B.  Winds,  Temperature,  Pressure  Altitude,  D-Values,  Ballistic  Winds  - 


HEIGHT (PRESS) 

DIR 

SPEED 

TEMP 

PRESS  ALT 

D-VALUE 

BALLISTIC  WINDS 

(KFT  AGL)  (MB) 

(DEG) 

(KNOTS) 

(DEG  C) 

(KFT  MSL) 

(FT) 

(DEG) 

(KNOTS) 

SFC(  959) 

93 

10 

20 

1.61 

+  0198 

93 

10 

1  (  926) 

101 

36 

22 

2.47 

+  0245 

96 

18 

2  (  894) 

119 

25 

25 

3.42 

+  0294 

108 

18 

3  (  863) 

129 

21 

24 

4.36 

+  0355 

108 

18 

Fig.  9  Example  of  a  JAAWIN  winds  forecast  file. 


The  program  will  automatically  recognize  what  type  of  file  was  chosen  and  how  to  read  it  properly.  The  wind 
data  will  then  be  loaded  into  the  appropriate  winds  matrix.  Auxiliary  data  such  as  the  name  and  type  of  the  file, 
when  and  where  the  winds  data  were  taken,  and  altitude  range  will  also  be  read  and  displayed  in  the  corresponding 
windows. 


B.  Drop  Zone  and  Map  Panels 

The  Drop  Zone  panel  allows  the  user  to  choose  one  of  the  drop  zones  stored  in  an  ASCII  file.  The  script  behind 
this  panel  loads  latitude,  longitude  and  mean-sea-level  altitude  of  a  chosen  drop  zone.  The  corresponding  map 
appears  on  the  left.  By  default  the  center  of  the  chosen  DZ  becomes  an  intended  impact  point  (shown  in  the 
corresponding  windows).  It  can  be  changed  by  selecting  a  new  point  on  the  map  or  typing  new  values  in  the 
windows. 

The  map  can  be  any  graphical  image  scaled  into  a  rectangle  with  respect  to  the  Universal  Transverse  Mercator 
(UTM)  coordinates.  When  the  map  loads,  the  axes  of  this  map  are  also  assigned  so  that  they  can  be  manipulated 
later.  For  instance,  the  following  set  of  commands  loads  the  jpg  image  of  some  map: 

A  =  imread( 'MAP_UTM_WGS84 ' , ' jpg' ) ; 

MAPimageR  =  A(;,;,l:3);  MAPimage ( : , ; , 1 : 3)  =  MAPimageR (end: -1 ; 1, ; , 1 : 3) ; 
clear  MAPimageR  A 

Mlat=[60.17  123.51] ; 

Mlon=[53.65  63.93]; 

image (Mlon,  Mlat,  MAPimage) ;  axis ( [Mlon  Mlat] ) 
set  (gca,  ' YDir normal ') ;  axis  equal 

The  user  can  zoom  the  map  in  and  out,  move  it,  measure  distances  between  selected  points,  get  the  coordinates 
of  any  point  selected  with  the  cursor,  and  read  the  coordinates  of  any  point  in  one  of  the  chosen  formats  (UTM  on 
the  map,  or  Lat/Lon  in  the  windows).  Two  functions  support  conversion  of  UTM  coordinates  into  Lat/Lon  and  back. 
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C.  Test  System  Panel 


The  Test  System  panel  allows  the  user  to  choose  one  of  the  ADSs  stored  in  an  Excel  file  (see  Fig.  10).  This 
system  can  have  one,  two,  or  three  stages.  The  corresponding  windows  for  unused  stages  are  shadowed.  The  only 
ADS  data  needed  for  fan  computation  are  the  system’s  descent  rate,  glide  ratio  and  altitude  to  deploy  (altitude  loss 
for  the  first  stage).  The  user  is  allowed  to  change  any  of  these  three  parameters  in  the  corresponding  windows  or 
create  an  entirely  new  system.  This  data  along  with  winds  aloft  and  release  point  information  is  used  to  compute 


safety  fans. 


Fig.  10  Examples  of  ADS  input  data. 


D.  Release  Point  and  Units  Panels 

The  idea  behind  the  Release  Point  panel  is  that  based  on  the  system  being  tested  and  available  winds  data,  the 
user  enters  various  release  points  and  computes  safety  fans.  If  the  safety  requirements  are  not  met,  the  user  moves 
the  release  point  and  a  new  safety  fan  can  be  generated.  The  release  point  can  be  moved  by  clicking  on  a  new 
location  in  the  map,  in  which  case  the  new  coordinates  immediately  appear  in  the  corresponding  panel  windows,  or 
by  entering  a  new  latitude  and  longitude  in  the  panel  windows.  Altitude,  airspeed,  and  direction  at  the  release  point 
are  provided  by  the  user.  When  a  satisfactory  release  point  has  been  chosen,  its  coordinates  will  be  provided  to  the 
crew  for  execution. 

As  an  example,  the  following  commands  are  activated  upon  clicking  the  “Move  Release”  button  in  the  GUI  and 
choosing  to  select  a  new  point  by  clicking  on  the  map.  The  current  release  point,  if  any  exists,  is  deleted  from  the 
map.  The  ginput  ( )  function  is  interactive  and  allows  the  user  to  click  anywhere  in  the  map  to  select  a  point.  The 
new  release  point  is  then  plotted  and  the  corresponding  Lat/Lon  calculated: 

%%  Deleting  current  release  point  if  any 

h  =  f indob j  (' MarkerFaceColor k -and MarkerSize 12 ) ; 
h2  =  findobj ( 'MarkerFaceColor r -and' , 'MarkerSize 6) ; 
delete (h);  delete (h2); 

%%  Getting  and  plotting  chosen  release  point 
[x, y]  =  ginput (1) ; 
rp  =  plot (x, y, ' k* ' ) ; 

set (rp,  'MarkerFaceColor  '  ,  '  k ' ,  'MarkerSize ' , 12 ,  ' LineWidth ' , 2 ) ; 
rp2  =  plot (x, y, ' ro ' ) ; 

set (rp2 , 'MarkerFaceColor ' , ' r ' , ' MarkerSize ', 6) ; 

%%  Determining  lat/lon 
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easting  =  x*10^3  +  350000; 
northing  =  y*10^3  +  1800000; 
zone  =  '  1  IS  '  ; 

[RPlatDeg, RPlonDeg] =utm2wgs (easting, northing, zone) ; 

%%  Accounting  for  sign  of  longitude  (using  West  in  GUI  vs.  negative) 

RPlonDeg  =  -RPlonDeg; 

Since  the  ground  crew  and  the  deployment  crew  can  speak  different  technical  languages  (units),  the  Units  panel 
allows  the  user  to  change  units  from  the  International  System  of  Units  (SI)  to  English  and  back.  In  addition,  Lat/Lon 
coordinates  can  be  visualized  throughout  all  panels  in  one  of  the  following  three  formats: 

•  Deg  -  degrees. decimal  degrees 

•  DM  -  degrees,  minutes. decimal  minutes 

•  DMS  -  degrees,  minutes,  seconds. decimal  seconds 

If  using  the  DM  or  Deg  formats,  the  unused  windows  are  shadowed.  For  example,  the  chosen  format  in  Fig. 5  is  DM, 
and  therefore  the  first  windows  in  all  Lat/Lon  readings  are  shadowed. 

E.  Operations 

Other  operational  functions  in  the  GUI  include  timing,  fan  computation,  loading  and  saving  data,  and 
commenting.  The  “Start  Clock”  button  asks  the  user  for  the  computer’s  time  zone  setting  and  starts  a  continuously 
updating  clock  in  Coordinated  Universal  Time  (UTC).  When  all  the  necessary  data  for  fan  computation  has  been 
selected  and  entered,  the  user  can  compute  the  current  safety  fan.  Corresponding  graphics  appear  on  the  map.  The 
time  of  the  most  recent  calculation  appears  under  the  running  clock. 

The  GUI  can  be  saved  with  the  Save  Data  button.  Once  a  GUI  has  been  saved,  it  can  be  loaded  with  the  Load 
Data  button.  When  reloaded,  the  GUI  is  still  interactive.  This  means  that  though  it  can  be  used  as  a  simple  figure 
capturing  one  state  in  time,  it  can  also  be  used  as  the  starting  point  for  continuing  operations.  If  using  a  loaded  GUI, 
the  clock  must  be  restarted  but  all  other  data  are  available. 

The  user  can  enter  comments  on  a  particular  configuration  in  the  Comments  panel.  These  comments  will  be 
saved  with  the  GUI.  The  Stop  GUI  button  stops  all  activity  in  the  figure  and  is  used  just  prior  to  closing  the 
application. 


VI,  Conclusions 

Computation  of  the  safety  fans  for  testing  of  aero  delivery  systems  involves  simple  and  yet  cumbersome 
computations.  For  uncontrolled  systems  the  precomputed  ballistic  winds  can  reduce  some  of  the  computational 
burden.  However,  for  effective  operations  with  any  system,  especially  the  ones  with  several  stages,  some  kind  of 
automation  is  required.  This  paper  presented  an  example  of  such  an  automated  tool  which  performs  several  safety 
checks  in  convenient  graphical  form.  The  mathematical  foundation  for  safety  fan  computations  was  also  rigorously 
developed.  The  graphical  user  interface  is  ready  and  will  be  used  soon  in  a  field  environment.  Feedback  will  be 
provided  to  the  developers  for  further  improvement  and  modifications. 
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